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Abstract. 

We present results of a study of the formation and evolution of the 
dark matter (DM) halos in a COBE-normalized spatially flat ACDM 
model (n = 1 - Q A = 0.3; h = 0.7). The dynamics of 256 3 DM par- 
ticles is followed numerically in a box of 60/i _1 Mpc with the dynamic 
range of 32, 000 in spatial resolution. The high resolution of the simu- 
lation allows us to examine evolution of both isolated and satellite halos 
in a representative volume. We discuss the new halo finding algorithm 
designed to identify halos in high-density environments, present results 
on the evolution of velocity function of DM halos and compare it with 
the Press-Schechter function, discuss the evolution of power spectrum of 
matter and halo distributions, and mass evolution of halos. 

The velocity function of halos at z = compares well with the pre- 
diction of the Press-Schechter approximation, but for circular velocities 
in the range 100 - 200 km/s simulations predict ~ 1.3 time more halos 
(mostly in clusters or groups). In real space the power spectra of halos and 
DM are very different (halos are anti-biased). Both spectra do not have 
simple power-law shape. In redshift space the spectra are close to a power 
law with 7 = —2.1 in the range of wave numbers k = 0.2 — 5/iMpc . 
The power spectra of halo distribution evolves only mildly for z = — 3. 
The mass evolution of isolated virialized objects determined from the 
simulation is in good agreement with predictions of the extended Press- 
Schechter models. However, satellite halos evolve very different: for some 
of them the mass decreases with time, which happens if the halos fall into 
clusters or groups. We discuss the dependence of the correlation function 
of halo populations on their environment and merging history. 
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1. Introduction 



A variety of observations on a wide range of spatial scales (i.e. rotation curves 
of galaxies, galaxy dynamics and gravitational lensing in galaxy clusters, the 
large-scale velocity flows, etc.) indicate the existence of dark matter (DM) in the 
Universe. The total amount of DM is not yet known, but it is generally believed 
that cold dark matter dominates the mass in the Universe and significantly 
affects both the process of galaxy formation and the large-scale distribution of 
galaxies. The gravitational domination of DM on the scale of galactic virial 
radius implies that collisionless simulations can be used to study formation of 
the DM component of galaxies. 

It is well known that cosmological scenarios with cold dark matter alone 
cannot explain the structure formation both on small and very large scales. Re- 
cently, scenarios with a non-zero vacuum energy, quantified by the cosmological 
constant A, have been proved to be very successful in describing most of the 
observational data at both low and high redshifts. In this contribution, we dis- 
cuss the evolution of DM halos in a spatially flat cosmological model dominated 
by the cosmological constant and cold dark matter, ACDM. A simulated box of 
60 h^ 1 Mpc size allows to construct and study statistically representative halo 
catalogs (about 10,000 halos are identified at z = 0). The novel feature of this 
analysis is study of both isolated and satellite halos, where halos are dubbed 
satellites if they are located within the virial radius of a larger host system 
(massive galaxies or galaxy groups and clusters). This allows us to study the 
differences between isolated and satellite halos; here we discuss the differences 
in their mass evolution and spatial distribution. 

2. Numerical simulation 

We simulate the evolution of 256 3 cold dark matter particles in a spatially flat 
model with a cosmological constant (ACDM; = 1 — f^A = 0.3; <7 8 = 1.0; 
Hq = 70 km/s/Mpc). The model is normalized in accord with the four year 
COBE DMR observations (Bunn & White 1997) and observed abundance of 
galaxy clusters (Viana & Liddle 1996). The age of the universe in this model is 
w 13.5 Gyrs. 

In order to study the statistical properties of halos in a cosmological envi- 
ronment, the simulation box should be sufficiently large. On the other hand, to 
assure that halo survival in clusters, the force resolution should be < 1 — 3/i _1 kpc 
and the mass resolution should be < 1O 9 /i _1 M0 (Moore et al. 1996, Klypin et al. 
1998, hereafter KGKK). The simulations were done using the Adaptive Refine- 
ment Tree (ART) code (Kravtsov, Klypin & Khokhlov 1997). The code used a 
512 3 homogeneous grid on the lowest level of resolution and six levels of refine- 
ment, each successive refinement level doubling the resolution. The sixth refine- 
ment level corresponds to a formal dynamical range of 32,000 in high density 
regions. Thus we can reach in a 60 h^ 1 Mpc box the force resolution of ~ 2/i _1 
kpc. The simulation has a mass resolution (particle mass of 1.1 x 10 9 h~ 1 M.Q) 
sufficient to identify galaxy-size halos (M>3 x 10 10 h~ 1 M.Q). 

The ART code integrates the equations of motion in comoving coordinates. 
However, its refinement strategy is designed to effectively preserve the initial 
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physical resolution of the simulation. In order to prevent degradation of force 
resolution in physical coordinates, the dynamic range between the start (zi = 30) 
and the end of the simulation should increase by (1 + Zi): For our simulations 
it should have at least the dynamical range 512 x (1 + Zj) = 15, 872. This 
is accomplished with the prompt successive refinements in high-density regions 
during the simulations. The peak resolution is reached by creating a hierarchy 
with six levels of refinement. The spatial refinement is accompanied by the 
similar refinement of the integration time step. The integration step of Aao = 
0.0015 corresponds to 645 time steps on the lowest-resolution of the uniform grid 
and to (effective) 41,280 on the deepest refinement level . In physical units, the 
latter step corresponds to 2.3 x 10 5 years. 

3. Challenges in halo identification 

Identification of halos in dense environments and reconstruction of their evo- 
lution is a challenge. Any halo finding algorithm has to deal with difficult 
"decision-making" situations, existing also in the real Universe. The most typi- 
cal difficulties arise if a small satellite is bound to a larger galaxy (like the LMC 
and the Milky Way or the M51 system) or in cases when many small gravita- 
tionally self-bound halos are moving within a large region of virial overdensity 
(galaxy clusters and groups). 

Assuming that the satellite is self-gravitationally bound, we would have to 
include the mass of the satellite in the mass of the host system. By doing so, 
we count mass of the satellite twice: when we identify the satellite and when we 
identify the host system. This may seem unreasonable, but if we do not include 
the satellite, then the mass of the large galaxy is underestimated. For example, 
the binding energy of a particle at the distance of the satellite will be wrong. The 
problems arise also when we try to assign particles to different halos in the effort 
to find halo masses. This is very difficult to do for particles moving between 
halos. Even if a particle at some moment is bound to one of the halos, it is 
not guaranteed that it belongs to the halo. The gravitational potential changes 
with time, and the particle may end up falling onto another halo. These effects 
are especially frequent in high-density crowded regions typical of galaxy groups 
and clusters. This situation actually was observed very often in simulated halos 
when we compared particle contents of halos at different redshifts. It represents 
an additional difficulty in tracing the mass history of such halos. 

It is also difficult to estimate the actual mass of a halo orbiting inside a 
cluster-size object. The formal virial radius of such a satellite halo is simply 
equal to the cluster's virial radius. We have chosen to define the outer tidal 
radius of the satellite halos as a scale at which their density profile starts to 
flatten. At small distances from the center of the satellite the density steeply 
declines, but then it flattens out and may even increase. This means that we 
reached the outer border of the satellite. The corresponding mass within tidal 
radius is then our approximation to the halo mass. With these considerations in 
mind, the halo finding algorithm must allow halos to overlap and DM particles 
to belong to more than one halo. 
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4. Halo finding algorithm 



The most widely used halo-finding algorithms, the friends-of-friends (FOF) and 
the spherical overdensity, both discard "halos inside halos" , i.e., satellite halos lo- 
cated within the virial radius of larger halos. The distribution of halos identified 
in this way, cannot be compared easily to the distribution of galaxies, because 
the latter are found within larger systems. In order to cure this, we have devel- 
oped two related algorithms, which we called the hierarchical friends-of-friends 
(HFOF) and the bound density maxima (BDM) algorithms (KGKK). 

Since the algorithms work on a snapshot of the particle distribution, they 
tend to identify also small fake "halos" consisting of only a few unbound particles, 
clumped together by chance at the analyzed moment. We deal with this problem 
by both checking whether the identified clump is gravitationally bound and by 
following the merging history of halos. A halo that does not have a progenitor 
at the previous moment is discarded, if the particles from which this halo has 
been formed belong at the previous moment already to other halos. In this case 
the halo is assumed to be a fake one. If it has formed from single particles or 
small objects below the threshold of halo detection, it is assumed that a new 
halo has been formed. For other halos we find the direct progenitor, i.e. a halo 
at a previous moment that contains the maximum number of particles of this 
halo. We use the chain of progenitors identified in this way to reconstruct the 
mass evolution of the given halo back in time, down to the epoch of its first 
detection in the simulation. 

The HFOF and BDM algorithms are complementary. Both of them find 
essentially the same halos. Therefore, we believe that each of them is a stable 
algorithm which finds in a given dark matter distribution all dark matter halos 
above a given mass threshold. The advantage of the HFOF algorithm is that it 
can handle halos of arbitrary, not only spherically symmetric, shape. The ad- 
vantage of the BDM algorithm is that it describes better the physical properties 
of the halos due to the fact that it separates background unbound particles from 
the particles gravitationally bound to the halo. 

5. Velocity function 

Interacting halos exchange mass and lose mass. The total mass of a halo depends 
on its radius which, as was pointed out above, is difficult to define. We try to 
avoid this problem by assigning not only a mass to a halo, but finding also 
its maximum "circular velocity" (>/ GM/R), v c i rc . This is the quantity which 
is more meaningful observationally. Numerically, v c i rc can be measured more 
easily and more accurately then mass. 

The output of the halo finding algorithm depends primarily on the assumed 
mass threshold. With the threshold of lO lo /i _1 M (10 particles) the algorithm 
identifies ~ 17000 halos, whereas a mass threshold of 3 x lO lo /i~ 1 M (28 par- 
ticles) results in identification of more than 9000 halos. There is a weak de- 
pendence on the assumed maximum halo radius (2% decrease if changing the 
maximum radius from 100/i _1 kpc to 150/i _1 kpc). 

For any study one needs to have a complete halo sample that is not affected 
by the numerical details of halo finding procedure. In order to test the complete- 



4 




100 1000 100 1000 



Figure 1. Differential (left panel) and integral (right panel) velocity 
functions of halos at z = 0. The thick lines show the velocity function 
of halos, while thin lines represent the Press-Schechter prediction. 

ness of the halo samples, we have constructed the differential velocity functions 
at z = for different mass thresholds and maximum radii. We find that the 
halo samples do not depend on the numerical parameters of the halo finder for 
halos with i> c j rc >100 km/s. Particularly, the differential velocity function for 
Vdrc > 100 km/s is robust and does not depend on the assumptions. There 
is a substantial scatter for v c i rc > 500 km/s where we have less than 10 halos 
per logarithmic velocity bin. In fact, increasing the maximum halo radius one 
(slightly) increases also the maximum circular velocity, if it is not yet reached. 
This is the case for the most massive halos. We reanalyzed these halos and found 
that in all cases the maximum circular velocity changed by less than 10 %. 

In Figure [[], we compare the velocity function of halos in our simulation 
with the Press-Schechter (PS) prediction. The latter has been been computed 
by converting virial mass in the PS mass function to v c i rc , assuming the NFW 
density profile of the halos. There is an excess of small halos (v c i rc < 200 km/s) 
over the PS prediction. This is expected because we include satellite halos in 
our catalog, while PS model predicts mass function of the isolated objects. For 
higher circular velocities we find a slight deficit of massive halos (20% - 30 %) in 
comparison with the numbers predicted by the PS mass function. In particular, 
for the largest halos, the over-prediction may be influenced by the assumption 
of NFW density profile which is a rather poor description for some of the large 
halos with substantial central sub-structure. The prediction of the simulation 
is also influenced by cosmic variance in our single realization; we may have a 
deficit of massive halos in the given realization just by chance. 

6. Halo evolution 

In the following we will discuss first the evolution of halo clustering and then 
the evolution of the halos. We consider three catalogs of halos with v c i rc > 120 
km/s at z = 0, 1, 3 (4804, 8018, 8396 halos respectively). These samples are 
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Figure 2. Evolution of the integral velocity function of halos. The 
solid, dot-dashed, and dashed curves correspond to z = 0, z = 1, and 
z = 3, respectively. 

are complete for z = 0, 1: the differential velocity function does not show any 
turnoff (see Figure at these velocities. We are missing part of the halos with 
Vdrc < 130 km/s at z = 3. To study the evolution of the halos we construct the 
mass evolution history of all halos which have been identified in the simulation 
at z = 0. 

6.1. Evolution of halo clustering 

In Figures |3| and ||, we show the evolution of the power spectrum of these halos in 
comparison with the power spectrum of the DM particles. We have calculated 
the power spectrum on a 512 3 grid; the spectra are shown up to half of the 
Nyquist frequency h^ ox x 128. The shot noise level S = V/N for the halo samples 
is 45, 27, and 26, respectively. The shot noise power derived from a random 
distribution is constant up to about a quarter of the Nyquist frequency kb ox X 64 
and decrease than by about 25%. We expect an error of the same order in the 
halo power spectrum at k^ ox x 128. 

The halo power spectrum evolves much slower than that of the dark matter. 
The halos are clearly biased at z = 3, are essentially unbiased at z = 1, and are 
anti-biased at z = 0. A detailed analysis of the evolution of halo and matter 
power spectra and bias will be presented elsewhere (Kravtsov and Klypin 1998, 
in preparation). The results obtained from the power spectrum analysis are 
in qualitative agreement with results on the bias evolution as derived from the 
correlation function analysis by Cohn et al. (1998) and other researchers. The 
scale dependence of the bias, however, is quite different. 

In the right panels of Figures || and ||, we show the power spectrum of halos 
and dark matter in redshift space. The redshift space power spectra are always 
much steeper then in real space. At z = 0, the real space power spectrum of dark 
matter shows a clear excess over the linear power spectrum due to nonlinear 
clustering. The corresponding dark matter power spectrum in redshift space 
follows almost exactly the linear power spectrum of dark matter. 

In real space the power spectra of halos and DM are very different. Both 
spectra do not have a simple power-law shape. The real-space power spectrum 
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Figure 3. Evolution of real-space (left) and redshift-space (right) 
power spectrum of the DM particles. The solid, dot-dashed, and dashed 
curves correspond to z = 0, z = 1, and z = 3, respectively. The dotted 
line left shows the linear power spectrum at z = 0. 




Figure 4. Evolution of real-space (left) and redshift-space (right) 
power spectrum of halos with circular velocity v c i rc > 120 km/s. The 
solid, dot-dashed, and dashed curves correspond to z = 0, z = 1, and 
z = 3, respectively. The dotted lines show the DM power spectrum at 
z = 0. 
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of halo distribution evolves only mildly between z = 3 and z = 0. Similarly to 
the real space case, the redshift-space power spectrum of halos shows almost no 
evolution during this time. In the range of wave numbers k = 0.2 — 5h Mpc -1 
the redshift space power spectrum of halos is close to a power law with 7 = —2.1. 
A power spectrum slope of ~ —2.1 had been measured for the combined SSRS2 
+ CfA2 galaxy sample on scales <30/i _1 Mpc (da Costa et al. 1994). 

6.2. Progenitors of halos 

In order to study the evolution of individual halos, we need to construct a 
complete evolution tree for each of the halos in the z = catalog. We proceed 
as follows: We have selected 19 epochs (z = 0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 
0.8, 0.9, 1.0, 1.5, 2., 2.5, 3., 5., 10., 15.). For every epoch, we identify a progenitor 
of a z = halo. The procedure of progenitor identification is based on the 
comparison of lists of particles belonging to the halos at different moments both 
back and forward in time. As was mentioned above, the halo finder algorithm 
allows halos to overlap, or, in other words allows particles to belong to more than 
one halo. The visual inspection of a large number of constructed evolution trees 
showed that this forward-backward algorithm of tracing halo histories identifies 
the correct "ancestor-descendant" relationships rather accurately, with obvious 
ancestor-descendant misidentifications in < 2% of the cases. 

Using the procedure described above, we are able to address the question 
of the halo detection epoch, which we define as an epoch at which the halo has 
been identified for the first time. The epoch of the first identification depends, 
of course, on the mass threshold assumed by the halo finder. The following 
results are for the lowest possible mass 3 x lO lo /i- 1 M at z = 0. For this mass 
threshold the halo catalog is dominated by halos of mass J; 10 11 /i _1 Mq, which 
fall below the identification threshold quickly as we trace their mass evolution 
back in time. The distribution of detection time for these halos is thus likely to 
reflect selection function determined by the threshold rather than any kind of 
physical "formation epoch" distribution. In the left panel of Figure ||, we show 
the distribution of detection epochs for the halos identified with mass threshold 
of 3 x lO lo /i -1 M . The number of progenitors at a given redshift exponentially 
decreases with redshift. The figure shows that the first halo in this catalog 
should have been detected as early as z = 10. Indeed, at z = 10 we find one 
halo of mass 1.95 x 10 11 M Q (184 bound particles shown as filled circles in the 
right panel of Figure |5|) . The overdensity of the halo is 400 which means that it 
is a virialized object. We show also the surrounding, not yet bound, particles of 
the halo (open circles). The overdensity in the box of 400/i _1 kpc size centered 
on this object is ~ 50. The figure also shows a flattened structure of the particle 
distribution: the halo forms inside of a small pancake. 

6.3. Halo mass evolution and environment 

The mass of an object found by the HFOF algorithm at virial overdensity can 
be defined as the sum of linked particle masses. In this case we do not only 
find galaxy size halos but also all group size and cluster size halos. For all of 
the HFOF objects we identify the main progenitors at all epochs down to the 
halo detection time. To study the mass evolution due to merging and accretion 
we have divided these objects into three mass bins at z = 0. The average 
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Figure 5. Left: Number of halos in a 60h Mpc simulation box 
whose progenitors existed at redshift z. Only halos with mass larger 
than 3 x lO lo /i _1 M have been taken into account. Dotted line: cx 
exp(-z/zi), z\ = 1.1. Right: A virialized halo (M 1.95 x 10 n M o ) 
at z = 10; filled circles show particles bound to the halo, the open 
circles show unbound particles. 



mass evolution of halos in these bins normalized to the mass at z = is shown 
in Figure |(| left. The dotted line is for average halo mass of 1.2 x 1O 13 M (12 
halos), the dashed for 1.0 x 1O 12 M (17 halos), the dash-dotted for 1.1 x 1O U M 
(265 halos). The solid lines show the predictions of semi-analytical model (kindly 
provided by Claudio Firmani) for 10 13 , 10 12 , and 1O 11 M , from the bottom to 
the top. We find a good agreement with the semi-analytical predictions (Lacey &; 
Cole 1993) for the evolution of the FOF selected objects, which are per definition 
isolated. 

Unfortunately, there is no simple and straightforward way to assign a mass 
for all halos identified in the simulation. Unlike the isolated halos identified by 
HFOF at virial over density, the satellite halos, although surviving, are subject 
to tidal stripping which reduces their mass. They are limited therefore by tidal, 
rather then virial, radius. To assign masses to the halos we proceed as follows. 
The isolated halos are assigned the mass inside the virial radius or radius of 
lOO/i -1 kpc, whichever is smaller. The satellite halos are assigned the total 
mass of gravitationally bound particles within their tidal radius (or, again, within 
lOO/i^ 1 kpc, whichever is smaller). The tidal radius is determined as the radius 
at which the density profile of a halo flattens (stops decreasing). 

We now construct the complete mass evolution histories for the set of all 
halos with the masses assigned as described above. We have divided these halos 
into five groups with masses M > 1O 13 M , 1O 13 M > M > 5 x 1O 12 M , 5 x 
1O 12 M > M > 1O 12 M , 1O 12 M > M > 5 x 1O U M , and 5 x 10 n M o > M . 
We defined a subset of 3674 halos, mass of which increases (with allowance for 
small statistical fluctuations) at all epochs. As before, the mass of these objects 
is normalized to their final mass at z = 0. The mass evolution of these halos is 
shown in Figure ||, right (solid lines). The solid lines are for average masses of 
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Figure 6. Left: Mass evolution of isolated objects identified by the 
HFOF algorithm at virial overdensity. The thin lines are for different 
average masses of these objects (see text). The solid lines show the 
predictions of the extended Press-Schechter approximation for different 
masses (see text). Right: Mass evolution of all halos. The solid lines are 
for different average masses (see text), the dot-dashed and the dashed 
lines show the mass evolution of a subset of halos which loose mass due 
to the tidal interaction. The radius of the halos was restricted to be 
not more than 100 h~ l kpc. 

(from the bottom to the top) 1.2 x 1O 13 M (14 halos), 6.6 x 10 12 M Q (34 halos), 
1.9 x 10 12 M Q (442 halos), 7.0 x 10 n M Q (534 halos), and 2.4 x 1O U M (2650 
halos). The overall evolution is similar to the mass evolution of isolated halos 
described above Figure |6[ left). Note, however, that while the mass evolution 
tracks are curved in Figure g, left, the mass evolution of the sample that includes 
satellites can be better represented by the straight lines in these log- log plots. 
This difference is due to the different halo selection procedure and to the different 
assignment of mass to the selected halos. 

In the two lowest mass ranges we also find an additional subset of 2650 
halos, whose masses decrease after z = 1: The dot-dashed (average mass of 
6.9 x 1O 11 M ) and the dashed (average mass of 2.0 x 1O 11 M ) lines in Figure |6] 
(right) show the mass evolution of this subset of halos which loose mass due to 
the tidal stripping in groups and clusters. 

Their mass increases at high redshifts, reaches a maximum and decreases 
thereafter. Contrary to the main halo population, the mass of which always 
increases {merging halos) the mass of these stripped halos grows first due to 
accretion of surrounding material and of smaller halos. At some point, however, 
these halos start to loose mass due to either tidal stripping or interaction with 
other nearby halos. In the left panel of Figure fj], we show the distribution of 
these halo subsets at z = 1 and z = 0. At z = 1 the "stripped" halos are 
distributed similarly to the rest of the halos. At z = 0, however, they are 
clustered much more strongly than the overall halo population. 
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Figure 7. Left: A slice of 15/1" 1 Mpc thickness: Merging (left) and 
stripped (right) halos at z = (top row), the same halos at z = 1 
(bottom row). Right: The correlation function of halos with circular 
velocities v c i rc > 120 km/s at z = 0. The solid line corresponds to halos 
the mass of which always increases (solid lines in Figure |6|, right), the 
dashed line corresponds to halos which loose mass during evolution 
(dashed and dashed-dotted lines in Figure ^ right). The correlation 
function of dark matter particles is shown by the dotted line. 

In the right panel of Figure 0, we show the correlation function for the two 
subsets of halos: always increasing mass and decreasing mass at z < 1. The 
correlation functions of the former has a lower amplitude and is not as steep as 
the correlation function of the latter. Note that the CF of the halos with the 
ever-increasing mass is anti-biased at scales < lOh^ 1 Mpc, the CF of the halos 
that loose mass is actually positively biased. This reflects the fact that the 
loosing mass halos are found within massive systems such as massive galaxies, 
groups, and clusters, and are therefore strongly clustered. 

One might speculate that this difference in the correlation functions may 
serve as a possible explanation for the color segregation of the correlation am- 
plitude that has been recently observed (Carlberg et al. 1998). In fact, one 
could expect that the galaxies hosting halos which undergo different mass evolu- 
tion also show different properties, and colors in particular. Further studies are 
necessary to test whether this simple picture can really explain the observations. 
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